import sys
import numpy as np
import matplotlib as MPL
import matplotlib.pyplot as pyplot
from scipy.misc import comb
from scipy.optimize import fsolve
from scipy import integrate as I

j = np.complex(0,1)

class Bernstein(object):
    
    Plots = []
    hold = []
    
    def __init__(self, points=[]):
        self._P = np.array(points)
        self._n = np.size(self._P, 0) - 1# the zero implies that the data must
                                          # be provided in the rows of 'points'
        self.beenThru = 0
        
    def s(self, t, i=0):
        ''' 
        Provides reparameterization of t from i to i+1 when self is the i-th 
        (0-based indexing) curve in a path; necessary when integrating along a
        piecewise continuous path; the default i=0 results in s = t, i.e., the
        standard parameterization of a single Bezier.
        '''
        return t - i

    def b(self, t, j=-1, i=0): # b(t;j)
        if self._n > 0:
            n = self._n 
            if j in range(n+1): # This means only b_j(t;i) is required
                return np.sum([comb(n,k) * np.power(1-self.s(t,i),n-k) * \
                                            np.power(self.s(t,i),k) \
                           * self._P[k][j]#.reshape(-1,1) \
                           for k in np.arange(n+1)],0)
            else:
                return np.sum([comb(n,k) * np.power(1-self.s(t,i),n-k) * \
                                            np.power(self.s(t,i),k) \
                           * self._P[k].reshape(-1,1) \
                           for k in np.arange(n+1)],0)
    
##    def bminusx(self, t, x): # used by tGivenCoord to find root t 
##                               # for guesses x when self._n > 2
##        b = self.b(self.s(t))
##        bmx = b - x
##        return bmx
##        return self.b(self.s(t)) - x

    def bprime(self, t, i=-1):
        if self._n > 0:
            n = self._n
            if i in range(n+1): # This means only b'_i(t) is requiredk
                s = self._P[n][i]*np.power(self.s(t),n-1) - \
                    self._P[0][i]*np.power(1-self.s(t),n-1)
                s += np.sum([self._P[k][i] * comb(n-1,k-1) * \
                        (np.power(1-self.s(t),n-k) * np.power(self.s(t),k-1)
                        - (n-k)*np.power(1-self.s(t),n-k-1) * \
                        np.power(self.s(t),k)/k)
                     for k in np.linspace(1,n-1,n-1)], 0)
                return n * s
            else:
                s = self._P[n].reshape(-1,1) * np.power(self.s(t),n-1) - \
                    self._P[0].reshape(-1,1) * np.power(1-self.s(t),n-1)
                s += np.sum([self._P[k].reshape(-1,1) * comb(n-1,k-1) * \
                      (np.power(1-self.s(t),n-k) * np.power(self.s(t),k-1) \
                        - (n/k-1) * np.power(1-self.s(t),n-k-1) * \
                        np.power(self.s(t),k))
                     for k in np.linspace(1,n-1,n-1)], 0)
                return n * s

    def B(self, range=np.linspace(0,1,101)): # b(range of t)
        return self.b(range)
    
    def Bprime(self, range=np.linspace(0,1,101)): # b'(range of t)
        return self.bprime(range)
    
    def SB(self, a, b, i=-1): # Integral of B
        B = lambda t,j: self.b(t, j)
        S = lambda j: I.quad(B, a, b, j)[0]
        if i == -1:
            SB = np.array([S(j) for j in range(len(self._P[0]))])
        else:
            SB = S(i)
        return SB
    
    def length(self):
        ds = lambda t: abs(np.dot(self.bprime(t).T, (1,j)))
        s = I.quad(ds, 0, 1)
        return s[0]
        
    def tGivenCoord(self, x, i): # i--0 or 1--specifies which coordinate x is
        if self._n==1:
            t = (x - self._P[0][i]) / (self._P[1][i] - self._P[0][i])

        elif self._n==2:
            disc = self._P[1][i]**2 - self._P[0][i]*self._P[2][i] + x * \
                        (self._P[0][i] - 2*self._P[1][i] + self._P[2][i])
            if disc < 0:
                t = np.NAN
            else:
                disc = np.sqrt(disc)
                A = self._P[0][i] - self._P[1][i]
                B = self._P[0][i] - 2*self._P[1][i] + self._P[2][i]
                t = ((A + np.array((disc, -disc))) / B)
                t = t[(t >= 0) & (t <= 1)]
                
        elif self._n==3:
            onethird = 1.0/3
            root3j = np.sqrt(3)*j
            u = np.array((1, (root3j-1)/2, (-root3j-1)/2))
            x0, x1, x2, x3 = self._P[:,i]
            a = x3 - 3*x2 + 3*x1 - x0
            if a != 0:
                b = 3 * (x2 - 2*x1 + x0)
                c = 3 * (x1 - x0)
                d = x0 - x
                del0 = np.complex(b*b - 3*a*c)
                del1 = np.complex(2*b*b*b - 9*a*b*c + 27*a*a*d)
                C = np.power((del1 + np.sqrt(del1*del1-4*del0*del0*del0))/2,
                              onethird)
                t = -(b + C*u + del0/C/u)/a/3
                t = t.real[np.isclose(t.imag,0)]
                t = t[np.abs(t-0.5) <= 0.5]
                if len(t) == 1:
                    print "tout: ", t
                else:
                    cb = self.b(t)
                    
                    print "b(t): ", self.b(t)
                self.beenThru = self.beenThru + 1
                return t
        else:
            pass
##            guesses = np.linspace(0.05, 0.95, self._n + 1)
##            ts = np.round(fsolve(self.bminusx, [guesses, guesses], 
##                                            (np.tile(x, guesses.size),)), 8)
##            print "ts: ", ts
##            tsy = np.round(fsolve(self.bminusx, guesses, (x[1],)), 8)
##            print "tsy: ", tsy
##            tmp = np.tile(x, guesses.size)
##            ts = np.round(fsolve(self.bminusx, guesses.T, (tmp,)), 8)
##            print "ts: ", ts
##            t = np.intersect1d(tsx,tsy) # just gets rid of repeats
##            print "t: ", t
##            print "self.b(ts): ", self.b(ts)

    def Plot(self, range=np.linspace(0,1,101), newFig=True, hold=False):
        if newFig | (not Bernstein.Plots):
            newFig = pyplot.figure()
            Bernstein.Plots.append(newFig)
            Bernstein.hold.append(hold)
        curFig = Bernstein.Plots[-1]
        curhold = Bernstein.hold[-1] & hold
        if curhold:
            if not hold:
                pyplot.gca().hold()
        x, y = self.B(range)
        pyplot.plot(x, y, axes=pyplot.gca())
        curFig.show()
        return newFig
    
                    
class Path(object):
    def __init__(self, Bs=[]):
        self._Bs = Bs

    @property
    def N(self): return len(self._Bs)

    def whereConnected(self):
        whereConnected = np.zeros(self.N, dtype=bool)
        if self.N:
            for n in range(self.N):
                prev = n - 1
                whereConnected[n] = (self._Bs[prev]._P[-1] == \
                                                    self._Bs[n]._P[0]).all() 
        return whereConnected
        
    def isClosed(self): # Note: assumes Bs are "well-ordered"
        return self.whereConnected().all()
    
    def isInside(self, P, iterTol=1e-9, resPrec=9): # uses complex contour
                                                      # integration 
        if self.isClosed():
            while (iterTol < 1e-5) and (resPrec > 5):
                S = 0
                for path in self._Bs:
                    z = lambda t: np.imag(np.dot(path.bprime(t).T, [1,j]) \
                                 / np.dot((path.b(self.s(t)).T - P), [1,j]))
                    absErr, newAbsErr = (1,0.5)
                    nSubintvls = 50
                    while (nSubintvls < 1e3) and \
                           (newAbsErr > iterTol) and \
                           (absErr != newAbsErr):
                        absErr = newAbsErr
                        tmp = I.quad(z, 0, 1, limit=nSubintvls)
                        newAbsErr = abs(tmp[1])
                        nSubintvls *= 2 
                    S += tmp[0]
##                print S % (2*np.pi)
                test = S % (2*np.pi)
                if np.isclose(test, 0) or np.isclose(test, 2*np.pi):
##                    print np.round(S % (2*np.pi), resPrec)
                    if -np.log10(iterTol) == resPrec:
                        iterTol *= 10
##                        print 'iterTol: ', iterTol
                    else:
                        resPrec -= 1
##                        print 'resPrec: ', resPrec
                else:
                    break
            if (np.isclose(iterTol,1e-5)) and (resPrec==5):
                S = np.NAN
            else:
                S = np.round(S, resPrec) != 0
        else:
            S = 0
##        print S
##        print
##        return S
    
    def isOn(self, p, tol=1e-6): # returns True if, to within tol, p is on 
                                   # self
        pass
        
    def isNorth(self, p): # True if p is north of self
        pass
        
    def isEast(self, p): # True if p is east of self
        pass
        
        